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Abstract 

The influence of random site dilution on the critical properties of the two-dimensional 
Ising model on a square lattice was explored by Monte Carlo simulations with the 
Wang-Landau sampling. The lattice linear size was L = 20 — 120 and the concen- 
tration of diluted sites q = 0.1,0.2,0.3. Its pure version displays a second-order 
phase transition with a vanishing specific heat critical exponent a, thus, the Harris 
criterion is inconclusive, in that disorder is a relevant or irrelevant perturbation for 
the critical behavior of the pure system. The main effort was focused on the specific 
heat and magnetic susceptibility. We have also looked at the probability distribu- 
tion of susceptibility, pseudocritical temperatures and specific heat for assessing 
self-averaging. The study was carried out in appropriate restricted but dominant 
energy subspaces. By applying the finite-size scaling analysis, the correlation length 
exponent v was found to be greater than one, whereas the ratio of the critical ex- 
ponents (a/u) is negative and (7/V) retains its pure Ising model value supporting 
weak universality. 
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1 Introduction 



In the last few years, many theoretical, numerical and experimental investiga- 
tions have appeared for the study of the influence of disorder (usually quenched 
site or bond dilution) on statistical systems in the presence or absence of an 
external magnetic field, since an experimental sample is not free of disorder, 
in general. Such systems are randomly dilute uniaxial antiferromagnets, e.g., 
Fe q Zni_ q F 2 , Mn q Zni_ q F 2 , Fei_ q Al q , obtained by mixing uniaxial antiferro- 
magnets with non-magnetic materials. These systems are modelled by using 
pure systems models modified accordingly. The most popular models for pure 
systems are those of Ising and Potts; in the current case, the former one shall 
be used. 

The Ising model, a favorite for physicists, is used as the prototype for phase 
transitions, critical phenomena, biological and econophysics, owing to the fact 
that its two-dimensional version on a square lattice was solved analytically 
by Onsager, [1]; it has also become a standard model for testing scaling and 
universality hypotheses. The Ising model as well as that of Potts, Heisen- 
berg, Baxter- Wu, etc., represent collective phenomena which are difficult to 
be solved exactly, except in some cases, consequently approximate methods 
have been developed to cope with, such as mean field theory, numerical meth- 
ods, perturbation theory, scaling, renormalization group, Monte Carlo etc. The 
effect of randomness on the critical behavior and magnetic phase diagrams of 
classical random spin systems has attracted much interest in recent years. 
Randomness is encountered in the form of vacancies, variable bonds, impu- 
rities, random fields, etc. Its influence on the critical behavior and magnetic 
properties is a long-standing and still unsettled problem in statistical physics. 
Similar methods, as in the pure systems, have been utilized for the study of 
these systems. In these studies, an important question which had arisen was 
the extent to which randomness influences the critical behavior and magnetic 
properties. The first remarkable criterion for the influence of randomness on 
a critical system was proposed by Harris, [2], according to which randomness 
changes the critical behavior if the specific heat critical exponent a of the 
pure system is positive, a > 0. In this case a new critical point with con- 
ventional power law scaling and new exponents emerges. However, the pure 
two-dimensional Ising model (2D IM) constitutes a marginal situation since 
a — 0. This exception makes 2D IM of particular interest attracting much 
attention. However, this effort led to contradicting results, increasing confu- 
sion and revealing the inhibit hidden complexity, implying the need for more 
subtle approaches to tackle it. 

The 2D IM consists of an array of N fixed points lying on the sites of a two- 
dimensional lattice of linear size L such that N = L 2 . Each lattice site is occu- 
pied by a magnetic atom characterized by the spin variable Si, i — 1, 2, . . . , N, 
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with Si taking on the values ±1. One can also consider a modified version 
of this model wherein some of the lattice sites, chosen randomly, are either 
vacant or occupied by non-magnetic atoms, (e.g., Al atoms [3]); both cases of 
randomness are treated equivalently. This version of Ising model is called two- 
dimensional randomly site-diluted Ising model (2D RSDIM). For this system, 
the critical exponents associated with the random fixed point have been esti- 
mated for the dilution-type disorder by theoretical and numerical approaches, 
leading to questionable conclusions, [4] . However, these approaches have shed 
light on estimates of the critical temperature and nature of the phase transi- 
tion. In addition, the phase diagram is strongly influenced, among other fac- 
tors, by dilution, [5] . The concentration of vacancies (or nonmagnetic particles) 
is denoted by q (dilution), while that of occupied sites (magnetic particles) by 
p (purity), q + p = 1. The vacancies are considered to be quenched and un- 
correlated, since one can also encounter systems where the vacancy locations 
are correlated, [6,7]. 

The process to be followed here for tackling the RSDIM is that of Monte 
Carlo. The MC approach has been proved to be a powerful tool to study dif- 
ficult problems such as random spin systems. In some cases, the simulation 
method suffers from problems of slow dynamics, thus, new algorithms have 
been proposed to overcome such difficulties. Wang-Landau (WL)[8] and en- 
tropic sampling [9] are examples of such efforts. The critical properties concern 
always an infinitely bulk system, since the phase transitions appear in such 
a system; however, from MC simulations the critical behavior is extracted 
from results obtained on a finite-size system by means of finite-size scaling 
(FSS). This process, since its inception, has been evolved as a very powerful 
tool for extracting properties of an infinite system near a phase transition, 
although it is, as yet, not fully completed causing, sometimes, ambiguities 
about its results. The major goal of the finite-size method is to identify the 
set of critical exponents that, together with other universal parameters, char- 
acterizes the universality class. As these exponents offer the most direct test 
of universality, their precise calculation is of great importance. However, the 
experimental devices are finite, consequently, the exponents cannot be mea- 
sured with infinite precision causing, occasionally, controversies to distinguish 
the universality class a specific system belongs, [10]. For more on the theory 
of finite-size scaling see, e.g., Barber [11], Privman [12] and Binder [13]. 

For the 2D RSDIM, the nature of the possible phase transition as well as 
the universality class is not completely understood. The value and even the 
sign of the specific heat exponent a is still not known. Because of logarith- 
mic divergence of specific heat of the pure system (a = 0), Harris criterion 
is inconclusive, hence, a great deal of effort has been dedicated to elucidating 
the properties of the 2D RSDIM. Currently, it seems that two scenarios have 
prevailed, which, however, are mutually exclusive. According to the first one 
(logarithmic- corrections scenario), the critical exponents are unaffected by dis- 
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order, apart from possible logarithmic corrections (strong-universality) , while 
the second, predicts critical exponents varying continuously with disorder, but 
the exponents' ratios (7/V), (P/v) remain the same as in the pure case, weak 
universality, [4]. 

Kim and Patrascioiu [14] studied the 2D RSDIM on a square lattice with pe- 
riodic boundary conditions and dilution-site concentration (dilution probabil- 
ity) Q = 1/9, 1/4, 1/3 and lattice linear size L up to 600, using MC simulation; 
their conclusions are that the specific heat does not diverge for q = 1/4, 1/3 
while for q = 1/9 it seems to increase as e — > 0, e = (T — T C )/T C ,T C critical 
temperature. The magnetic susceptibility x an< l correlation length £ fit the 
pure power law, the value of the respective exponent 7 and v increases with 
q while r\ = 2 — 7/// remains the same as in the pure system. Queiroz and 
Stinchcombe [15] using transfer-matrix-scaling technique, Mazzeo and Kiihn 
[16] following the same technique with the equilibrium ensemble approach to 
disordered systems, came to the same conclusions. Newman and Riedel [17], 
using renormalization group on weakly diluted systems confirmed the exis- 
tence of a new stable fixed point with new exponents. Heuer [18] focused on 
the exponents (7/^) and ((3/v) for the 2D RSDIM on a square lattice with 
72 < L < 250 and < q < 0.4; according to his estimations, {(3/v) does not 
change notably with dilution within the errors and, practically, it assumes its 
pure system value. The exponent (7/^) shows a strong dilution-dependence in 
the temperature range 10~ 2 < (T — T c )/T c < 1, but near T c it asymptotically 
approaches the pure system value (7/4) independently of dilution. 

On the other hand, Shchur and Vasilyev [4] analyzing the MC data for the 
magnetic susceptibility critical amplitudes (r, V) above and below the critical 
temperature, respectively, for dilution q < 0.25 and lattice linear size up to 
L = 256, concluded that the ratio (T/V) seems to remain constant for the 
dilute-site concentrations considered and equal to its pure system value. This 
implies that the 2D RSDIM belongs to the same universality class as the 
pure one. Dotsenko and Dotsenko [19], Shalaev [20], Shankar [21], Ludwig 
[22], using field theoretical calculations, showed that randomness is irrelevant 
(critical exponents are unchanged) and only logarithmic corrections might 
appear in the case of weak dilution; they found for the correlation length, 
magnetization and susceptibility, respectively, 
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where A is a smooth function of q, a = and C a constant. Ballesteros 
et al [23], by performing MC simulations in conjunction with finite-size scal- 
ing (FSS) for p = 1,8/9,3/4,2/3, although observed small deviations of the 
values of the critical exponents from those of the pure system on varying 
the concentration, they considered it as a transient effect, since this can be 
lifted if a pure Ising value for the exponents is combined with logarithmic 
corrections. Selke et al [24], using MC techniques for lattices with linear size 
8 < L < 256 and spin concentration 0.1 < p < 1, concluded that impurities 
lead the specific heat to diverge as C ~ ln(lnL) on approaching the critical 
temperature. Tomita and Okabe [25] performed MC simulations on the same 
system on a square lattice using the probability-changing cluster (PCC) al- 
gorithm confirming that randomness is irrelevant and its influence is evident 
through logarithmic corrections. 

Allowing for the previous arguments, we remark that in spite of much effort 
devoted to the investigation of RSDIM in the critical region, the question of 
the dependence of critical exponents on randomness is still open. In this paper, 
we shall examine the critical properties of the 2D RSDIM using an alternative 
approach based on the WL algorithm. In this analysis, the major goal will 
be to estimate the critical temperature and exponents of the 2D RSDIM for 
various values of the spin concentrations to assess whether it belongs to the 
Ising universality class or not by studying the finite-size behavior of the specific 
heat and magnetic susceptibility. 

The paper is organized as follows. In the next section, after the introduction of 
the model, we shall discuss the approach based on the WL algorithm and an 
efficient implementation by conveniently restricting the simulation in the dom- 
inant energy subspace. In section 3 we apply the FSS to the Ising model under 
consideration and we close with the conclusions and discussions in section 4. 

2 Numerical approach of the RSDIM 

We consider the Hamiltonian of the two-dimensional site-diluted Ising model, 
in the absence of any external field, 

H = -J CiCjSiSj, Si = ±1, (5) 

<ij> 

where J > is the interaction constant, ferromagnetic interactions. The coef- 
ficients Ci's, called occupation variables, are quenched, uncorrelated random 
variables chosen to be equal to 1 with probability p, when the i-site is occupied 
by a magnetic atom and with probability q = 1 — p otherwise; that is, we 
have the probability distribution -P(q) = pS(ci — 1) + q5(ci). The summation 
extends over all nearest-neighbor pairs of the square lattice, of linear size L 
with periodic boundary conditions. 
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The data was generated by extensive MC calculations using WL sampling 
method to estimate the density of states (DOS) g(E) [8]. WL sampling per- 
forms a random walk with an acceptance ratio P(Ei — > Ej) = min{l, [g(Ej)/g(Ej)]}, 
Ei and Ej are the energies before and after the transition, respectively, aim- 
ing at sampling a flat histogram in energy. The WL algorithm overcomes the 
difficulties, such as critical slowing down and long relaxation times in systems 
with complex energy landscape, appearing in other MC processes. This algo- 
rithm has been applied to a wide range of systems spanning from the Ising 
model [26,27], random field Ising model [28], 3D conserved-order-parameter 
Ising model [29], Potts [30], to glassy systems [31], polymers [32], DNA [33] 
and the Baxter- Wu model [34]. The DOS g(E) is not constant during the 
random walk, it is modified according to the rule g(E) — > {g{E) ■ /); the mod- 
ification factor / varies as fj + \ = \J~fj, j is the order of iteration. In the current 
investigation, the WL algorithm performed 26 iterations and the initial value 
for / was f — e. 

Having an accurate estimation of g{E), the non- normalized canonical distri- 
bution can be constructed, P(E,T) = g(E)e~ f3E ; subsequently, the partition 
function can be calculated through the expression Z(T) = J2e g(E)e~ l3E , from 
which most of the thermodynamic observables can be estimated. This kind 
of MC simulation constitutes a major improvement towards speeding calcu- 
lations since it avoids multiple runs (one for each temperature) needed by 
the majority of MC algorithms to describe the temperature dependence of 
thermodynamic quantities over a significant temperature range; in a WL al- 
gorithm temperature is not needed to be specified a priori. In the present 
case, the WL algorithm was implemented on lattices with 20 < L < 120 and 
the density of states was stored as a function of the energy. The dilution q 
can vary from 0.0 (p = 1) to the percolation threshold q EERC = 0.407255 
(pPERC = 0.592745(2)), [35]. 

The quantities, upon which we shall rely for studying the critical behav- 
ior of the RSDIM, are the specific heat (C = ([< E 2 >] - [< E >] 2 )/{NT 2 )) 
and magnetic susceptibility (x = ([< M 2 >] — [< M >] 2 )/(NT)), per parti- 
cle. The MC data is generated by choosing a realization of the dilution for a 
specific value of q and various values of L. This procedure is then repeated for 
several other realizations for a specific value for the lattice linear size L. For 
each realization, we find the maximum value of the specific heat (C*(q, L)) and 
susceptibility L)), recording simultaneously the values of the respective 

pseudocritical temperatures, T^(q,L) and T*(q,L), respectively, forming in 
each case, a sequence of "pseudocritical temperatures" converging to the crit- 
ical temperature T c (q) of the infinite system as L — > oo. Among the different 
realizations large fluctuations are observed in the same set of the previous 
quantities. 

The presence of randomness is evident in the way averaging processes are 
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carried out for an observable X, which assumes a different value for each of 
the M random realizations of the disorder corresponding to the same value of 
dilution q. This implies that X behaves as a stochastic variable, whose mean 
value is estimated through a two-step process. First, the usual thermal average 
is performed for a specific realization of the dilution. After the completion of 
the M realizations of randomness, the disorder average is carried out over the 
M realizations and is denoted by the brackets [ ]. 

For the 2D IM the specific heat and magnetic susceptibility in the thermo- 
dynamic limit diverge; in a finite lattice, this divergence is rounded off and 
manifests itself by a maximum exhibited by the above quantities. This maxi- 
mum increases gradually with L and ultimately tends to infinity as L — > oo. 
For the 2D RSDIM, in attempting to detect the maximum C*(q,L) of the 
specific heat and the respective temperature T*(q,L) (pseudocritical temper- 
ature), we considered two routes. Let C m (q,L) be the specific heat for a par- 
ticular realization m out of M realizations for a specific dilution q. In the 
first route, we estimated, initially, the maximum value C^(q, L) together with 
the respective pseudocritical temperature T5 m (g, L) for every realization m of 
the disorder; then, we considered the sample average of the individual specific 
heat maximum [C*{q,L)\ and pseudocritical temperature [T£(q,L)], for the 
M realizations, 

1 M I M 

[C*(q,L)\ = - £ C* m (q,L), [T*(q,L)] = - £ T*, m (q,L) (6) 

m=l m=l 



In the second route, following Rieger and Young, the sample summation for 
the specific heat for the totality of M realizations was considered [36] , 

■y m 

[C(q,L)] sum = -J2 C rn(q,L) (7) 

m=l 



In (7), the resulting specific heat curve is very complex with many local max- 
ima, reflecting the strong pseudocritical temperature fluctuations in the en- 
semble of random realizations. From these maxima we selected the absolute 
one, indicated by [C(q, L)]* um = max[C(q, L)] sum , occurring at the pseudocrit- 
ical temperature T£< sum (q, L). The same procedure was also followed for the 
magnetic susceptibility x\ the resulting quantities are denoted by [x*(q,L)], 
[T;(q,L)],[x(q,L)Y sum ,Tl sum (q,L). 

To investigate the critical behavior of the system, we performed extensive MC 
simulations to calculate the density of states g(E) for each value of dilution q 
and L through the WL algorithm. After estimating the density of states g{E), 
one can proceed to the calculation of the necessary thermodynamic quantities, 
such as the energy E, specific heat Cl(T), magnetization M and susceptibil- 
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ity Xl(T) for further use in order to identify the probable universality class. 
However, before proceeding to this end, we outline a new restricted method 
for speeding up the numerical calculations; this method is called " Critical 
Minimum Energy Subspace" (CrMES) technique, [27,29,34]. In this method, 
we concentrate our simulation only on the dominant energy subspaces. In the 
following lines we describe its implementation. Consider the specific heat per 
site for a lattice of linear size L at temperature T, 

{Emax 
Z- 1 £ E 2 exp[S(E)-(3E] - 

E ■ 

[Z- 1 ]T Eexp[S(E)-pE}) ] (8) 

E 

where the Boltzmann constant was set ks = 1, thus (3 = 1/T, d is the spatial 
dimension (d = 2) and Z the "partition function" of the system, 

Z = E ^exp[S(E)-(3E] (9) 

E 

The latter expression is the partition function in case g(E) is the exact DOS 
of the system and properly normalized, [37]. In practice, the DOS, resulting 
from WL simulations, is an approximate result whose accuracy depends on 
that of the simulation. In the expression (8), the calculation of the specific 
heat in the critical region can be speeded up by restricting the energy in- 
terval if we use the CrMES technique. Let E be the energy corresponding 
to the maximum term exp[S(E) — (3E\ of the partition function (9) for the 
temperature at hand. Because of the sharpness of the energy distribution, the 
energy interval (E m i n , E max ) in the summation (8) is replaced by a smaller one 
(E_,E + ) around E corresponding to a predefined accuracy r for the specific 
heat expressed as, | [C L (E-, E + )/C L (E min , E max )} - 1 |< r, where r = 1 • 1(T 6 
and E± = E ± A±, A± > 0. The induced errors are much smaller than the 
ones in determining the DOS, for more see [27]. The magnetic properties were 
obtained using the final stages of the WL algorithm, following our earlier 
proposal as in [38]. 

3 Finite-size scaling analysis. Results 

The FSS is based on the assumption that the free energy of a system of linear 
size L and in the absence of an external magnetic field scales as, 

E(L,e) = L^E (eL e ) (10) 
where tp — (2 — a) /v. The scaling of the correlation length £ = £ e^ suggests 
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that 9 = v 1 . The scaling function Fq(x) is universal, in that, it is independent 
of the lattice size. 




I 1 I 1 I 1 I 1 I 1 I 1 I 1 I 1 

1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 
InlnL, InL 



Fig. 1. Specific heat finite size approximations against InL (squares) and InlnL 
(circles) for q = 0.3. 

The quantity that specifies the scenario the RSDIM shall follow is, never- 
theless, the specific heat. According to the InlnL scenario the specific heat 
diverges as in (4), while according to weak universality and FSS, the specific 
heat obeys the power law, 

[C*(q,L)]=p 1 + qi L a ^ (11) 

In this expression the lack of the correction terms speeds up the convergence 
and yields more stable fits. 

The choice of the suitable asymptotic law for the specific heat has caused 
considerable concern, since one has to choose between Eqs. (4) and (11); thus, 
we have checked both. In the pure two-dimensional Ising model, q = 0, the 
specific heat against InL is a straight line, representing the forthcoming diver- 
gence, while against InlnL the respective curve bends upwards, see Fig. 6 in 
[16] as well as Ref. [39]. Plotting the specific heat data for q = 0.3 against InL 
and InlnL, respectively, see Fig. 1, we observed that the former curve (right- 
most in Fig. 1) bends downwards deviating from a straight line, whereas the 
latter curve (leftmost) seems to bend also downwards implying that the data 
does not follow the InlnL scenario. The same behavior was also observed for 
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Table 1 

The critical temperature and critical exponents resulting from the sample average 
of the specific heat (6) and the respective expression for susceptibility, for different 
values of spin dilution q. For each q-value, the first line corresponds to the specific 
heat data and the second to the susceptibility. 



ajv 



C c 



0.1 1.90668(0.00101) 
1.90271(0.00219) 

0.2 1.50675(0.00137) 
1.50180(0.00342) 

0.3 1.15280(0.00385) 
1.10051(0.00211) 



1 i^209 03699 
i - iozuy 0.03476 

i i ri oi 0.03813 
1.10 181 .03576 

1 189nQ°- 02647 

1 189SQ0-04776 
iiiozoy .04419 

1 9QKOQ0. 05739 
l.ZODZy .05251 

1 9Q71Q0.02949 
l.zo/ lO .02815 



-0.26406(0.02552) 2.38683(0.09096) 



-0.30892(0.01719) 1.24453(0.01400) 



-0.38223(0.03119) 0.71851(0.00659) 



1.74968(0.01369) 



1.74918(0.00983) 



1.74939(0.01643) 



the dilution q = 0.2, but this is more evident for the stronger dilution q = 0.3; 
for q — 0.1 it is not so evident because this is a crossover case. 

Mazzeo and Kiihn ([16]) checking the credibility of both scenarios, studied 
initially the lnlnL one; by fitting their data to this law, they concluded that 
it was difficult to accept definitely this scenario (4); instead, they focused on 
the power law (11) and estimated the a-exponent for various dilutions, test- 
ing the validity of the respective values by using the hyperscaling relation in 
combination with the //-exponent, estimated earlier. Their data resulted from 
calculations on the equilibrium ensemble approach together with numerical 
transfer matrix technique, phenomenological renormalization group scheme 
and conformal invariance on finite-width strips. 

In addition, we have also invoked the x 2 -test to discriminate between Eqs. (4) 
and (11) by estimating the ratio between the x 2 -test for the lnlnL scenario with 
the one for the power law suggested by (11): for q = 0.2 this ratio is 2 while 
for q = 0.3 it is 3, indicating an increasing tendency with q; hence, although 
the above description makes evident the need for data on much larger systems 
sizes, we shall adopt the power law (11) as it seems to be more suitable and 
reliable for the present case. This practice shall be followed in the sequel. 

Firstly, we fit the specific heat sample average maxima to the power law (11). 
Considering a specific value for q and extrapolating towards L — > oo the 
respective values of [C*(q, L)}, one can read off the asymptotic value pi = C°° 
as well as the critical exponents ratio {pt/v)\ these appear in Table 1 for 
the dilutions q = 0.1,0.2,0.3. We have also considered that the maximum 
values [C(q,L)]* sum of the specific heat resulting form sample summation, see 
(7), follow a similar power law as in (11), but with different symbols for the 
respective quantities, 
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Table 2 

The critical temperature and critical exponents resulting from the sample summa- 
tion for the specific heat (7) and the respective expression for susceptibility, for 
different values of spin dilution q. For each q-value, the first line corresponds to the 
specific heat data and the second to the susceptibility. 



a/u 



C z 



0.1 1.89557(0.00218) 
1.90329(0.00322) 

0.2 1.49638(0.00272) 
1.50281(0.00387) 

0.3 1.08592(0.00621) 
1.07301(0.00426) 



1.15335[j;S -0.26412(0.02986) 2.29955(0.0988) 
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-0.30877(0.01983) 1.18223(0.01361) 



-0.38221(0.03831) 0.53077(0.00704) 



1.74934(0.00694) 



1.74931(0.00963) 



1.74925(0.02756) 




Fig. 2. Specific heat fitting according to the FSS prediction (11). The specific heat 
corresponds to the sample average (6) and site dilution q = 0.1. For q = 0.2,0.3 the 
plots are similar. The coefficients Pi,P2,P$ correspond to p\, q\, a/vc, respectively, 
see Eq. (11). 

[C(q,L)}: um =p 1 + q 1 L^o (12) 

Fitting the respective values to (12) and extrapolating towards L — ► oo, the 
values for pi and (a/vc) are estimated and appear in Table 2 for the same 
q values. In both Tables, the asymptotic values of the specific heat form two 
decreasing sequences as q increases and it seems that they tend to zero at the 
percolation threshold (q c = 0.407255); the non- divergence is also evident by 
the levelling off of the specific heat data, indicative of approaching a saturation 
value, see Fig 2; the specific heat follows a similar behavior for the other 
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two values of q, as well. The main outcome of this fit is that the exponents' 
ratio (ac/v) appears to be negative, exhibiting a steady decreasing tendency 
approaching the percolation value ajv = —0.5, [40]. 

For any value of q, the respective critical temperature T c is not known a 
priori (as it happens to be for the random bond counterpart [41]) although 
there exists a formula T c (q) = T C (0)(1 — 1.565g), T c (0) is the pure system 
critical temperature [42], but of limited applicability since it is valid only for 
small impurity concentrations, thus the estimation of the T c for any q is of 
great importance. The correlation length exponent v is directly related to the 
specific heat exponent a through the hyperscaling relation a + dv — 2 (d 
system's dimensionality, d = 2) which acts as a constraint; the nondivergence 
of the specific heat (a < 0) is consistent with v > 1 for the current model. 
Assuming the FSS prediction, 



[Tc{q, L)] — T CjC + b\L~ x l v ° 
[T;(q,L)]=T ca + c 1 L- 1 ^ 



(13) 
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Fig. 3. Specific heat pseudocritical temperature fitting according to the first equation 
of the FSS prediction (13). The sequence corresponds to the sample average (6) and 
site dilution q = 0.1. For q = 0.2, 0.3 the plots are similar. The coefficients P\,P2, P3 
correspond to T C) c,bi, (l/vc), respectively, see Eq. (13). 



one can estimate the critical temperatures T c c and T c x resulting from the 
sample averages of the specific heat and susceptibility, respectively, as well as 
the corresponding exponent u, see Fig. 3. Also, similar scaling laws are con- 
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Fig. 4. Susceptibility fitting according to the first equation of the FSS prediction 
(15). The specific heat corresponds to the sample average and site dilution q = 0.1. 
For q = 0.2,0.3 the plots are similar. The coefficients Pi,P2,Ps correspond to 
ro, si, respectively, see Eq. (15). 

sidered for respective quantities for the sample summation but with different 
symbols for the involved quantities, see [36], 

T^ sum (q,L) = f c , c + b 1 L- 1 ^ 

Tl sum (q ) L) = f c , x + Z 1 L- 1 ^ (14) 

Fitting the respective pseudocritical temperatures for both routes of the spe- 
cific heat and susceptibility to (13) and (14), they yield the values of the 
critical temperature and v exponent appearing in Tables 1,2. In both Tables, 
the respective values for v vary continuously with dilution q and are greater 
than one in conformity with hyperscaling. It seems that the critical tempera- 
ture decreases to zero as the dilution q increases towards the percolation limit 
q c = 0.407255. For the low impurity concentration, q = 0.1, the deviations 
of the respective values of the critical temperature is small, as well as for the 
intermediate q = 0.2, see Tables 1,2. For q = 0.1, the average of the four values 
is (T c (q = 0.1)) = 1.9020625, which agrees with that in Heuer [18] (1.9004427), 
Tomita and Okabe [25] (1.9022), Shchur and Vasilyev [4] (1.9032), and con- 
sistent with that from the formula in Ref. [42], T c (q = 0.1) = 1.9141. For 
the intermediate impurity concentration q = 0.2, the average of the respective 
critical temperatures is (T c (q = 0.2)) = 1.5016825; it is in agreement with that 
in Heuer [18] (1.507873), Shchur and Vasilyev [4] (1.5103), while deviates from 
that of the formula in Ref. [42], T c (q = 0.2) = 1.558930. For the large impu- 
rity concentration q = 0.3, the first value of critical temperature in Table 1 
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(1.15280) shows significant deviation from the other three, so if we consider 
only the remaining three their mean value is {T c (q = 0.3)) = 1.08648, consis- 
tent with that in Heuer [18] (1.075140), and Tomita and Okabe [25] (1.0712), 
while disagrees with that in Ref. [42] T c (q = 0.3) = 1.10306. 

We have also considered the FSS for both susceptibility averages. The maxi- 
mum values [x*(q, L)] and [x(q, L)]* um for both routes follow the scaling laws, 

[ X %q,L)}=r + s,L^ 
lx(q, L)\: um = r + 8 X tfP* (15) 

In (15), we have considered that both background terms vanish (ro = tq = 0), 
since this gives more stable fits. Fitting the respective numerical data for both 
routes of susceptibility to (15), the ratio (7/V) was estimated and their values 
appear in Tables 1, 2 for the same dilutions, see Fig. 4. In both Tables, the 
ratio (7/V) retains its pure Ising model value (7/V = 7/4) independently of 
dilution, thus corroborating the weak universality hypothesis. 




Fig. 5. Probability distribution of the relative susceptibility vs relative susceptibility 
Xr {Xr=Xj/(x)i Xji (x) are the j-bin and global average susceptibility, respectively) 
for dilutions q = 0.1(a), 0.3(6) and lattice linear size L=30, 70, 90. The non-self-av- 
eraging behavior is evident from the persistence of the width of the plots as the 
lattice linear size L increases. 

An important issue arising in a disordered system is the notion of self-averaging, 
that is, to what extent properties of the system depend on the particular re- 
alization of the quenched random variables, implying that the distribution of 
an observable becomes sharper as L increases. In case their relative widths 
remain constant as L — > 00, then we say the system is non-self- averaging, 
[43]. In the pure Ising model the respective distributions transform into delta 
functions as L — > 00. To investigate the possibility of such a behavior in the 
current disordered model, we produce N s samples of an observable X (e.g., 
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specific heat, susceptibility) and form the respective probability density dis- 
tribution P{Xj/ (X)), where Xj and (X) are the value of the observable in 
the j'-bin and its global average, respectively. If we identify X with system's 
susceptibility x, then the resulting probability density P(x%/ (x)) as a function 
of (xi/ (x)) f° r ? — 0.1, 0.3 and L = 30, 70, 90 appears in Fig. 5. The respective 
curves of the relative susceptibility for a specific value of dilution collapse on 
each other irrespective of the value for L and they do not become sharper on 
increasing L. The general shape of the curves remains the same independently 
of the value for q. This behavior implies non-self-averaging. Similar behavior 
was also displayed by the respective plots of pseudocritical temperatures and 
heat capacity. 

4 Conclusions and discussions 

We have presented MC numerical data on the 2D RSDIM for three cases 
of dilution spanning a wide range of dilution, q = 0.1,0.2,0.3, i.e., weak, 
intermediate and strong, for lattices with linear size in the range [20,120], 
using FSS and following the Wang-Landau algorithm for the calculation of 
the density of states. The calculations dealt with the estimation of the critical 
temperature, the critical exponent v and ratios (ayV) and (7/V). The results 
indicate that as q increases v increases, while (7/V) remains constant and 
(07V) decreases. A consequence of the invariance of (7/V) is that the exponent 
77 = 2 — (7/V) is invariant, as well. Using the Rushbrooke equality a+2/3+7 = 2 
in conjunction with hyperscaling relation a + du = 2, we deduced that (/3/V) 
is also invariant, retaining its pure Ising value, (3/v = 1/8. Referring to the 
critical temperature, although there exist four different sequences, resulting 
from the various procedures for any dilution q, each one exhibits a decreasing 
tendency, tending to zero as q tends to the percolation limit q c = 0.407255. 

A notable feature of the current-model data is the asymptotic behavior of 
the specific heat as L increases; it tends to a saturation value for any q, in 
contradistinction to its counterpart in the pure model that diverges. The other 
critical quantity, susceptibility, still diverges as L — > 00. 

An important issue is that if impurity concentrations larger than q = 0.3 to- 
wards the percolation limit q c are considered, then the high dilution shall re- 
duce significantly the long-range correlations enhancing in this way the finite- 
size effects; to moderate this effect larger lattices are needed at the cost of 
larger computer-execution time. 

Our results and findings are supporting the argument that the 2D RSDIM 
does not belong to the same universality class as the pure 2D Ising model but 
to a new one. In conclusion, the present results favor the weak universality 
scenario. 
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